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Abstract 

This paper investigates the use of multiple directions of stratification as a variance reduction 
technique for Monte Carlo simulations of path-dependent options driven by Gaussian vectors. 
The precision of the method depends on the choice of the directions of stratification and the 
allocation rule within each strata. Several choices have been proposed but, even if they provide 
variance reduction, their implementation is computationally intensive and not applicable to re- 
alistic payoffs, in particular not to Asian options with barrier. Moreover, all these previously 
published methods employ orthogonal directions for multiple stratification. In this work we 
investigate the use of algorithms producing convenient directions, generally non-orthogonal, 
combining a lower computational cost with a comparable variance reduction. In addition, we 
study the accuracy of optimal allocation in terms of variance reduction compared to the Latin 
Hypercube Sampling. We consider the directions obtained by the Linear Transformation and 
the Principal Component Analysis. We introduce a new procedure based on the Linear Approx- 
imation of the explained variance of the payoff using the law of total variance. In addition, we 
exhibit a novel algorithm that permits to correctly generate normal vectors stratified along non- 
orthogonal directions. Finally, we illustrate the efficiency of these algorithms in the computation 
of the price of different path-dependent options with and without barriers in the Black-Scholes 
and in the Cox-IngersoU-Ross markets. 

Keywords. Monte Carlo methods, variance reduction, stratification methods. 

I Introduction 

The main purpose of Monte Carlo (MC) simulations is to compute integrals numerically. It is fre- 
quently the only alternative for solving problems in applied sciences and notably for financial ap- 
plications. The pricing of derivative contracts and value-at-risk calculations for risk-management 
purposes typically require numerical simulations. However, the MC method for high-dimensional 
problems is a demanding computational task and a considerable number of studies have been de- 
voted to increase its efficiency via variance reduction techniques. This paper investigates the 
use of multiple directions of stratification as a variance reduction technique for MC simulations 
of path-dependent options driven by high-dimensional Gaussian vectors. The precision of the 
method depends on the choice of the partitions of the space and the allocation of the number of 
samples within each strata. Usually, the strata are polyhedrons delimited by hyperplanes orthog- 
onal to a few direction vectors. Several choices have been proposed: Glasserman et al. [8] select 
the directions for the stratification of linear projections based on the quadratic approximation 
of the integrand or payoff function. In contrast, Etore et al. [4] find the directions by adaptive 
techniques. These two approaches provide a high variance reduction but their implementation 
can be computationally intensive and the former one cannot be applied to more realistic payoff 
functions such as Asian options with barrier at each time step. Moreover, these two methods 
suppose orthogonal directions for multiple stratification. In this work, we investigate the use of 
algorithms producing convenient directions, generally non-orthogonal, combining a lower com- 
putational cost with a variance reduction that is comparable to the above mentioned methods. 
In addition, we study the accuracy of optimal allocation, combined with the above stratification 
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techniques, in terms of variance reduction, compared to "fixed" allocation procedures such as 
Latin Hypercube Sampling (LHS). We consider the directions produced by the Linear Transfor- 
mation (LT) decomposition introduced by Imai and Tan jg] and the Principal Component Analysis 
(PC A). Moreover, we propose a new procedure based on the Linear Approximation (LA) of the 
"explained" variance of the payoff function by the use of the law of total variance. Notably, we 
design a novel algorithm that permits to correctly generate multivariate normal random vectors 
stratified along non-orthogonal directions. We illustrate the efficiency of the proposed algorithms 
and their combination for the computation of the price of different path-dependent options with 
and without barriers in the Black-Scholes (BS) and in the Cox-Ingersoll-Ross (CIR) models. In the 
former dynamics, it turns out that the LA and the LT approaches return the same first order direc- 
tion while this vector is almost parallel to the one obtained by the GHS technique even in the case 
of Asian options with a barrier at expiry. This justifies the application and the good performance 
of the LA (and LT) if the barrier is at each monitoring time. Consequently, the approaches return 
the same variance reduction and the LA (LT) is easier to implement and has a lower computational 
cost. We repeat our numerical investigation in the CIR framework where we find explicit solutions 
for the LT and LA directions. In order to find a further direction, we compute the first principal 
component of the sampled covariance matrix of the price process obtained by a MC estimation 
via a pilot run. In both BS and CIR dynamics, LT and LA return remarkable variance reduction 
with a low computational cost. We also show that in some setting the stratification along multiple 
directions can be more efficient than stratifying along a single one. In particular, the combination 
of the LA (LT) direction and a non-orthogonal direction, notably the first principal component, 
can even outperform the variance reduction of two orthogonal directions in the case of barrier 
options. Finally, as far as the allocation of the samples is concerned, in any case the LHS displays 
a considerable higher computational time and has always a lower variance reduction as compared 
to the use of a convenient direction of stratification with optimal allocation. 

The paper is organized as follows. Section [2] reviews the main ideas of stratification and 
the motivations of this study. Section [3] presents the new algorithm that permits the stratification 
along non-orthogonal directions. Section [4] discusses the use of convenient stratification directions 
and in particular, contains the presentation of the LT decomposition and the introduction of the 
LA procedure. In Section [5] we explain the financial applications and find the explicit solutions 
for the LA and the LT methods both for the BS and the CIR dynamics. In Section |6] the variance 
reductions and the computational costs of the proposed technique are illustrated by numerical 
experiments. Finally, Section [t] concludes the paper by summarizing the most important findings. 



2 Stratified Sampling and Linear Projections 

Stratified Sampling is a general variance reduction technique that consists of drawing the observa- 
tions from specific partitions of the sample space. More specifically, suppose we want to compute 
by MC simulations an expectation of the form E[g(Y)] where g : IR"^ ]R is a Borel function and 
Y is a R*^ -valued random vector with the assumption that E[^(Y)^] < 00. Consider a stratification 

variable X and let Ai, . .., be disjoint subsets of the real line for which P (^Uf=i G ^/c}) = 1- 
Then 

E[g(Y)] = ^ E[^(Y)|X e A,]P(X e A,) = ^ ]E[g(Y)|X e A,]p, (I) 
where pjt = 1P(X G ^k)i ^ = 1, . . . , -K. The stratified estimator with Ng draws is defined as: 

Lni-Ls{yki) = ^E'iLsiykj), (2) 

k=l ;=1 fc^i Hk j^i 

where n/c are the number allocations in the fc-th stratum and qi^ = nj^ / Ns is their fraction in the 
k-th stratum and Y/^j are independent draws from the conditional distribution of Y given X G Aj-. 

Its variance is given by I^f^j Pl'ni' where cr^- is the conditional variance of ^(Y) given X G Aj^. 

This estimator may be more efficient than the usual MC sample mean estimator of a random 
sample of size Ng. The potential higher efficiency of the former estimator critically depends on 
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the allocation rule and the choice of the partition of the sample space. The optimal allocation rule 
is the one that minimizes the variance of the stratified sampling estimator given the partition of 
the state space and the constraint Y^=\ = 1- It is given by: 

% = =K ■ (3) 

E)c=l VkO-k 

The probabilities p)t are known whereas generally the conditional variances are not known. They 
can be estimated in a pilot run and then used in a second stage to determine the stratified estima- 
tor. This is not the optimal procedure and more sophisticated techniques can be employed, see 
for example Etore and Jourdain [[5]. 

We focus our attention on MC simulation driven by high-dimensional Gaussian vectors that 
are of particular interest in financial applications. As such, we consider in the following only 
normal random variables. 



2.1 Stratifying Linear Projections: i-dimensional Setting 

We begin with a general description of stratifjdng a linear projection of a Gaussian random vector. 
Suppose Z is a d dimensional centered Gaussian random vector, Z ^ M{Q, S,z) and then consider 
the stratification variable X as the linear projection of Z over a fixed direction v G , X = v ■ Z. 
X is also Gaussian with variance v ■ 2^zv. This choice permits to partition the sample space WJ^ 
into strata defined by 

S,t,. = {xeM'', x-ve (4) 

Due to the Gaussian structure of the random variables we can generate Z stratified along the 
direction v in the following way. Consider a general Gaussian random vector Y = (Yj, Y2): 

vMv„v.,^^( (;;.), (I; ^-)) (5, 

and denote £ (Y1IY2 = x) the law of Yi given Y2 = x, it is possible to prove (see for instance 
Glasserman |7j) that 

£(Yi I Y2 = x) = A/- (jii + ^12^22^ (x - 112) , Ell - 2:12^22^^21) . (6) 

where we assume that E22 is invertible. Adapting the above result for Z given X = v ■ Z and 
Var[X] = V ■ 2^2 V = 1 we have 

£ (Z |X = X ) = i^^x, Ez - ^zvv^I^z \ ^ _^ rs^vx,Ez - Ezvv^Ez) . (7) 
\v-Lzv V ■ Ezv / v / 

If we consider Ez = /rf the above equation becomes: 

£(Z|X = x)=A/'(v;c,Jrf-v^v). (8) 

The conditional covariance matrix D = — v^v does not depend on x and since D is an orthog- 
onal projection matrix, we have DD^ = D. Due to this result, we do not need to compute the 
Cholesky (or a general square-root) matrix of D to sample from the conditional distribution of 
Z given X. These observations give an easy and simple algorithm to generate K samples of Z 
stratified along the direction v. 

Suppose now that Aj^ is the interval between the quantiles of order and of order of the 
standard normal distribution. We can sample from Z given Z ■ v G Aj. in the following steps: 

1. generate U ^U{[0,1])- 

2. Set y = ^ and X = ^~^{V), with O the inverse of the cumulative normal distribution. 

3. Generate TJ ^ M (0, independent on U. 

4. Set vX (I - vv^) Z'. 

We suggest to implement the last term in the last step as Z' — v(v ■ Z') which requires 0(d) 
operation rather than O(d^). 
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2.2 Stratifying Linear Projections: Multidimensional Setting 

We start with the case of orthogonal directions and consider a matrix V G '^ '^^ ,d' < d, whose 
columns are the direction vectors, such that V^V = I^t. Following the notation introduced above 
we have: 

X = V^Z (9) 

where now X is d' dimensional. Moreover, 
Consequently 

£ (Z |X = x) = (^^zV (v^^zV) X, Ez - Ez V (v^^zV) V^Sz) (11) 

where we assume that V^Ez V is invertible. In the case Hz = have 

£{Z\X = x) =J\f (^V (v^Vy^ x,Ia-V {V^Vy^ V'^^ . (12) 

Hence, if we adopt orthogonal directions = I^i the algorithm to stratify Z given X = V^Z is a 
simple multidimensional version of the algorithm illustrated before where now we should stratify 
the d' dimensional hypercube [0, l]'''. Suppose, for example, that we stratify the j-th coordinate 
of the hypercube, ; = 1, . . . , d', into Kj intervals of equal length so that we have a total number of 
Ki X ■ ■ ■ X K^i equiprobable strata. In this multidimensional setting we can sample from Z given 

X = y^Z e Af,, where = H/li ^ ( ^/ Tj )'^^ following steps: 

1 . generate U = ( !Ji . . . , JJ^/ ) with independent components each of law W ( [0, 1 ] ) . 

2. Set Vj = with ;■ e {1, . . . , and A:^ e {1, . . . , Kj}. 

3. SetX={X^...,Xa>),Xj = ^~HVj). 

4. Generate Z' --^ A/^ (0, 1^) independent of U. 

5. Set VX + (k - W^) Z'. 

We now investigate the possibility to stratify over different directions that can be non-orthogonal 
either. When the directions are not orthogonal the components of X are not independent since 
Var[X] = ^ I^, and the previous multidimensional algorithm cannot be adopted anymore. 

A first way yo approach this problem may be to assume X = Cx£ with e ~ Id') independent 
on Z and Cx G R'''^'^' such that Var[X] = CxC^, and use the following slight modification of the 
above algorithm. 

1. generate U = (!Ji . . . , 11^/) with independent components each of law W( [0, 1] ). 

2. Set Vj = with ;■ = {1, . . . , d'} and fcy = {1, . . . , Kj}. 

3. Seie={ei...,ed,),ej = ^-\Vj). 

4. Generate Z' ^ M (0, 1^) independent of U. 

5. Set V{Cl)-'e +{ld-V {V^V) y^)Z'. 

However, although mathematically correct, this algorithm stratifies the marginals of the random 
vector e that has independent components. This construction does not consider the fact that the 
marginals of X are not independent and the introduction of the dependence can affect this partial 
stratification in complicated ways (see Glasserman [7]). 
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3 Stratification along non-orthogonal directions 

In this section we show how to generate multivariate normal random vectors, Z ^ Af{0,I^), 
stratified along non-orthogonal directions. We prove the following proposition: 

Proposition i. Let Bj = {ei . . ., e^/} be a set of linearly independent vectors in M''', d' < d, such that 
||e,|| = 1, i = 1, . . . , d' , let B' = {i[ . . . , f^,} be the set of orthogonal vectors produced the Gram-Schmidt 

procedure: f' = e,- "nc'w^ '" • finally consider = {f/ = irffir,' = 1, the orthonormal 

W^mW ll*il 

version ofB' and let F be the d x d' matrix whose i-th column is i^. 

Suppose g ■.'K'^ such that E [g^(Z)] < +oo and consider two vectors in WJ^' , = {a'^ , .. ., a'^}, 
such that a~ < a'^ , V/ = 1, . . . , d' . We have 

E[g(Z)|fl- <Z-e/<fl+, /■ = !,...,/] 



E 



m=l 



T(m-l) 



L7„ o u 



a+ (tt('«-1) 



P (flj" < Z ■ e/ < f = 1, . . . , d') 



where 



O ( g„ ( U 



(13) 



(14) 



with the notation, U^^^ = (Ui, . . . , iiy) , / = 1 . . . , n anrf . . . , JJrf/ /./.rf. uniformly distributed random 
variables, all independent on Z; we assume Uq = and af = af. 

Remark i. The above result requires the computation of the joint probability P (ar < Z ■ e; < a^, f = 1, . . . , d') 
if/zere i/ze random variables Z ■ e,-, / = 1, . . . , d' are not independent; in contrast, this term is not necessary 
for the estimation o/E [g(Z)]. Indeed, suppose K strata, by conditioning we have: 



E [g(Z)] = 12 E [g{Z) |Z e k-th stratum] P (Z e fc-t/t stratum) , 



(15) 



f/zen plugging in the conditional expectation the result of equation ( [T3I tfte probabilities at the numerator 
and at the denominator simplify out. 

Proof. For simplicity we suppose d' = 2, the Gram-Schmidt procedure returns = fi = ei, 
= 62 — (ei ■ 62)61 and £2 = pf|y. It follows that 



E 



a-^ < Z ■ 61 < fl 
flj^ < Z ■ 62 < 



E 



a. < Z ■ fj < fl^ 



"2 ^(ere2)erZ 

11/^11 



< Z ■ f , < 



a]'-(ere2)erZ 

11/^11 



Based on the results of the Section |2]and the properties of the conditional expectation, the previous 
expression equals: 



CE 
where 



z + fa (0-1 (o^-) + u, (0(4) - - fi ■ z) ) ^,-^u,)<i 



2 (Lri)<f2-z<4(Ui) 



(16) 



C 



P 



flj < Z ■ 61 < 
fl^ < Z ■ 62 < fl^ 
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The expected value is then: 



E 







g (Z + fi (O-I (O(fl-) + «i (0(4) - 0(a-))) - . z)) l,-(„,)<f,z<,,. 
g ( Z + fi (o-i (O(fl-) + Ml (0(fl+) - <!>(«-))) - fi ■ Z) 



1 

E 



+ f2 ($-1 (0(fl2 (Wi)) + 112 (^(^("l)) - ^i'^li^l)))) - f2 ■ Zj 
X (0(g+(Mi)) -0(82 (Ml)))] ^^"1 

g( Z + fi (0-1 (O(fl-) + Hi ($«) - <!>(«-))) - fi ■ Z 



= E 

+ f2 (0-1 (0(fl2-(lfl)) + U2 (0(fl2+(LIi)) - 0(fl2 (Lfl)))) - f2 

X (0(4(111)) -<l>(S2("i)))]. 

Rearranging the terms in Z we get equation ([T3J for d' = 2. The result for d' direction is obtained 
iterating the steps above. □ 

4 Convenient Directions 

Given an allocation rule, the crucial point in the stratification of linear projections is the choice 
of the directions of stratification. Indeed, stratified sampling eliminates the sampling variability 
across strata without affecting the sampling variability within strata. Good directions are charac- 
terized by their higher capacity to dissect the state space into strata where the integrand function 
is nearly constant. In the following we describe the approaches that we adopt in order to find the 
directions of stratification. 

4.1 Principal Component Directions 

Suppose we want to find the srngled-factor approximation of a d-dimensional Gaussian random 
vector X ^ M{0, S) that maximizes the variance of v ■ X. This is equivalent to the following 
optimization problem: 

arg max v ■ Sv (17) 

l|v||=l 

Suppose Ai > ■ ■ ■ > Arf represent the eigenvalues of Z, in increasing order, and ei, . . . , their as- 
sociated eigenvectors, then the optimization above is solved by v* = ei an eigenvector associated 
to the largest eigenvalue Ai . 

As ei produces the linear combination ei ■ X that best captures the variability of the compo- 
nents of X. We may choose this vector as the first direction of stratification. In the case we would 
consider multiple stratification, we can iterate the optimization above. This means that we would 
consider ey,; = 1, . . . ,d, associated to the /-th eigenvalue, as the ;-th direction of stratification. In- 
deed, in the statistical literature, the linear combinations ey ■ X,;' = 1, . . . ,d, are called the principal 
components of X. The variance explained by the first k < d principal components is the ratio: 

EtiA, 
Ef=iA,- 

Finally, we note that this procedure based on the PGA only produces orthogonal directions. 

4.2 Law of Total Variance and GHS Directions 

In this section we illustrate the law of total variance and we briefly describe the strategy to select 
optimal directions illustrated in Glasserman et al. [8J. Given two random vectors Xi and X2 of 
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dimension and dj, respectively, and a function g : R^^i M, if E[^(X)^] < oo, the law of total 
variance reads as: 



Var[g(Xi)] =E[Var[g(Xi)|X2]]+Var[E[g(Xi)|X2]]. (18) 

Usually, in the context of linear model, the two terms are known as the "unexplained" and the 
"explained" components of the variance, respectively. In our case, Xi is a standard normal ran- 
dom vector Z and X2 = v ■ Z where v G R . It is well known that stratification eliminates the 
"explained" component of the variance up to terms with order o{l/Ns), where Ng is the total 
number of draws (see for instance Lemma 4.1 in Glasserman et al. [SJ). Hence, a good direction 
candidate is the one that maximizes the "explained" component of the variance or minimizes the 
"unexplained" part. 

Such an optimal direction is then the solution of the following optimization problem: 



arg mrn / Var g{Z) v ■ Z = x px{x)dx, (19) 
veR'',||v||=liE'' L J 



veR'',||v||^ 

where px is the density of X = v ■ Z. 

The approach proposed in Glasserman et al. JSl is to adopt directions that are optimal for the 
quadratic approximation of the logarithm of the integrand function. Glasserman et al. 18J consid- 
ered ^(z) = exp ^jz ■ Bz^ with B non-singular symmetric matrix whose eigenvalues Aj, . . . , A^; 
are all less than 1/2. Now number the eigenvalues and eigenvectors of the matrix B so that 



i-MJ - \1-X2J - \i-K, 

Glasserman et al. fE\ proved that the optimal direction v* is the eigenvector ei of the matrix B 
associated with the eigenvalue Aj. When one considers multiple stratification, the j-th optimal 
direction is the eigenvector associated with the eigenvalue Aj. Since the directions are the 
eigenvectors of the matrix B, the GHS approach only produces orthogonal directions. 

When the logarithm of the integrand function is not quadratic, one could evaluate its Hessian 
at the certain point. Glasserman et al. ||8| proposed to calculate the Hessian at a point used for an 
importance sampling procedure. This last operation might be really computationally expansive, 
in particular if d is large. It depends on a non-convex optimization procedure and cannot always 
be easily applied to realistic situations arising in finance. In addition, in financial applications, 
payoff functions (integrand functions) are far to be quadratic. In contrast, Etore et al. [4] found the 
directions by adaptive techniques that in some cases outperform the above approach. However, 
the numerical procedure still remains computationally intensive. These drawbacks motivate our 
study where our main purpose is to investigate convenient multiple stratification directions that 
provide comparable variance reductions with a notable advantage from the computational point 
of view. 



4.3 Linear Approximations 

In this section we describe a different approach, that we name Linear Approximation (LA), in 
order to find convenient directions for the stratification of linear projections. 

Suppose g & C^, this approach is based on a linear approximation of the function g that leads 
to an approximation of the "unexplained" component of the variance. Then, we can approximate 
the optimization problem fig) as: 



/ 



Vg(0) ■ Var 



Vg{0)px{x)dx, 



(21) 



where we also use the approximation Vg(E[Z Z ■ v = x]) ~ Vg(E[Z]), that is we evaluate the 

gradient at the expected value of Z (zero for each component) instead of its conditional one. The 
solution of the optimization problem pJ^ is given by the following proposition: 
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Proposition 2. The optimal direction v* of the optimization problem ( I21D is: 



Vg(0) 
|V^(0)| 



(22) 



Proof. Developing equation pl\ we get: 



V^(0) ■ Var 



X 



Vg{0)px{x)dx 



Vg{0)-{I-v^v)Vg{0)px{x)dx 



l|V^(0)f-V^(0)-v^vVg(0). 



(23) 



The minimization problem is equivalent to maximize the second term that can be written as 
(Vg(0) ■ v)^. The maximum of this dot product is attained when the two vectors are parallel. The 
optimal direction is then obtained by normalization. □ 

Multiple directions in the LA procedure can be produced calculating the gradient at different 
points. For example, we might iteratively consider Z2 = Vg (V^(0)) , . . . , Z^/ = \7g (V^(Z(;/_i)) 
in order to capture higher order components. We remark that the LA approach does provide 
non-orthogonal directions. 



4.4 Linear Transformations 

The LT procedure, proposed by Imai and Tan fgl, is originally conceived to enhance the accuracy 
of simulation techniques that employ low-discrepancy sequences also known as Quasi-Monte 
Carlo (QMC) methods. Indeed, given Z ^ J\f{0, 1^), the variance of the MC estimation of the 
expected value E[g(Z)] does not change if we replace Z by Ae where e ^ A/'(0, /j) and A is 
a d X d orthogonal matrix, AA^ = Z^, while the choice of A can deeply affect the accuracy of 
QMC simulations (see for instance Papageorgiou |i4j). The Imai and Tan's choice is such that 
A minimizes the effective dimension in the truncation sense defined in Caflisch et al. [3] of the 
integrand function. In our context, the columns of A will be chosen as the orthogonal directions 
of stratification. 

We briefly describe the LT algorithm. Consider a d dimensional normal random vector T ~ 
Af{}i; E), a vector w = {zui, .. ., Wd) e R'' and let /(T) = X]f=i ^iTj be a linear combination of T. 

Let C be such that Z = CC^ and assume e ^ J\f{0, I^j) with T = Ce. The LT approach considers 
C as C = C^'^ = C'-^A, with C*-^ the Cholesky decomposition of Then, in the linear case, we 
can define: 

g^ie) := /(C™Ae) = ^ a.e, + p-xv, (24) 
k=l 

where aj.- = C^^ • w = A.^ • B, fc = 1 . . .,d and B = (C'-^)^w while C.^ and A.^ are the k-th 
columns of the matrix C and A, respectively. In the linear case, setting 

A*i = ±p|, (25) 

with arbitrary remaining columns with the only constrain that AA^ = I^, leads to the following 
expression: 

g^{e) = p-w±Me^. (26) 

This is equivalent to reduce the effective dimension in the truncation sense to 1 and this means to 
maximize the variance of the first component ei . 

In a non-linear framework, we can use the LT construction, which relies on the first order 
Taylor expansion of g^: 

g^ie)^g^{e) + {:^-^Aei. (27) 
1=1 "'^^ 
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The approximated function is linear in the standard normal random vector Ae ^ Af{0, 1^) and we 
can rely on the considerations above. The first column of the matrix A* is then: 



A.i* = arg max ( ) (28) 



Since we have already maximized the variance contribution for (^-^^:^ j / in order to improve 

the method using adequate columns we might consider the expansion of g about d — 1 different 
points. More precisely Imai and Tan [9] propose to maximize: 

A.k* = arg max^ f ^4^) (^9) 



subject to ||A.k*|| = 1 and A.j* ■ A.k* = 0,; = 1, . . - l,fc < d. 

Although equation I23I 1 provides an easy solution at each step, the correct procedure requires 
that the column vector A.^* is orthogonal to all the previous (and future) columns. Imai and 
Tan [9] propose to choose e — e\ = E[e] = 0, q = (1,0, . . .,0), . . .e/t = (1, 1,1, . . .,0, . . .,0), 
where the A:-th point has k — 1 leading ones. Sabrno [16] illustrated an economic and convenient 
implementation of the LT algorithm by an iterative QR decomposition that we will use to find the 
directions of stratification. This method is computationally more expensive than the LA and it is 
not clear if it admits a solution when the sequence of expansion points is different from the one 
described above. 



5 Financial Applications 

In this section we illustrate how to calculate the convenient directions introduced above in the 
context of option pricing. We consider two price-dynamics: 

• BS d5mamics for M risky assets with constant volatilities: 

dSi{t) = rSi{t)dt + aiSi{t) d]Ni(t) , S,(0) = S,o, z = l,...,M, (30) 

S, [t) denotes the z-th asset price at time t, (t, represents the volatility of the i-\h asset return, r 
is the risk-free rate, and W (t) = (Wj (t) ,. . . , (0) ^ M-dimensional Brownian motion 
such that dW,(f)rfWjt(i) = Pikdt, i,k = 1,...,M. WhenM= 1 we simply denote S(0 = Si{t). 

• CIR d5mamics: 

dS{t)=a{ii-S{t))dt + a^Js{t)dW{t), S(0) = Sq, (31) 

with Sq, a, ji, cr positive constants. We impose the condition 2a}i > in order to ensure that 
S{t) remains positive. 

Applying the risk-neutral pricing formula (see Lamberton and Lapeyre [12]), the calculation of 
the price at time t of any European derivative contract with maturity date T boils down to the 
evaluation of an (discounted) expectation: 

fl(0 =exp(-r(T-0)]E[i/'| J-f], (32) 

the expectation is under the risk-neutral probability measure and is a generic J-'j-measurable 
variable that determines the payoff of the contract. 

We show how to derive the convenient directions of stratification for the following derivative 
contracts: 

1. discretely monitored Asian basket options: 

a(t) = exp(-r(T- f))E 



^ M N ^ 
EL^rjSi{tj) -Ks 







)■ 





Option on a Basket (33) 
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where x+ — max{x,0), ti < t2 - ■ ■ < tjq — T is a time grid, the coefficients Wij satisfy 
Y^i j Wij = 1 and Ks is the strike price. When N — 1 and M > the option is known as 
basket option while if M = 1 and N > it is simply known as Asian option. 



2. Asian option with knock-out barrier at expiry T: 



a{t) =exp(-r(T-f))E 



N 



N 




S(T)<B 



.Ft 



where B represents the value of the barrier. 
3. Asian option with knock-out barrier at each monitoring time: 



a{t) =exp(-r(r-i))E 



N 



E S (tj) - 1 ls(f^)<B,V;=l,...,N 



where B represents the value of the barrier. 



(34) 



(35) 



5.1 Linear Transformation in the Black-Scholes Market 

Suppose the BS dynamics with constant volatilities and a time grid ti < t2 - ■ ■ < t}^ — T, the 

elements of the autocorrelation matrix Eg of the Brownian motion are (Sb)^,, = min(ty, t„), j, n = 
1,...,N. Moreover, denote the a covariance matrix whose elements are {'La)^ — ^iPim'^m, 
i,m — 1, . . . , M, and consider Emn = Z,b ® where denotes the Kronecker product. Given 
e ~ N(0, Imn) and C^^ = C'^^A such that C'^^{C'^^Y = Emn and AA^ = Imn, the payoff of an 
Asian basket option can written as: 



MN 



MN 



= (^(e) - ^^s)"^ where g{e) '^n^i 



fc=i 



and 



= ln(w)fcifc2%(0))+ r 



(36) 



(37) 



where the indexes k\ and ^2 are fci = (fc — l)moduloM + I,k2 — [ik — 1)/MJ + 1, respectively 
and [xj denotes the greatest integer less than or equal to x. 

Since the Asian payoff function is not everywhere differentiable, the LT procedure is applied 
to its differentiable part g {or g — Ks). This is done also for the other barrier-style Asian options, 
hence we obtain the same directions of stratification for the three types of derivative contracts. 
Hereafter we detail the adopted procedure: 



1. Expand g up to the first order: 



NM /NM ( NM \ \ 

8{e) = g{e) + E ( E exp U + E Q^^H Cfj Ae, 



(38) 



2. For e = find the first column of the optimal matrix A: 



NM /NM \ 

gie) = g(0) + E ( E ^^P (Fi) Aei 



(39) 



Set ai = (ESf exp (f,) Cf) = E~5i {zlt exp (fO Cf^) ^mi and set u^) = (e^ 1 emN)T 

and B(i) = (C™)^u(i) then the first coliram is 



A*i = ± 



bW 

iB(i)ir 



(40) 
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3- The p-th optimal column is found considering the p-th expansion point of the strategy. This 
results in: 

NM /NM / p-i \ \ 

gie) = g{ip) + E ( E exp im + E C*A Cfj Aei (41) 

where C*j. = (C'-^A|);, k < p have been already found at the p — 1 previous steps and A.p* 
must be orthogonal to all the other columns. 

Also define u(P) = (exp {^fi^ + E^i Qc) exp (?^mn + ECi' ^l^^k) ) ^ b(p) = (C™)^u(P), 

then the solution is 

B(p) 

A* = ± — r^. (42) 

P ||B{P)|| 

We remark that at each time step all the columns must be orthogonalized (see Sabino [ [13} [16] ') 
5.2 Linear Transformation in the CIR Market 

We extend the procedure described in the previous section with the assumption of a CIR dynam- 
ics. Consider an equally spaced time-grid whose time step is denoted by At, the Euler scheme of 
the CIR dynamic is: 



Sj = Sj_i + (x[ii- Sj^i) At + a^Sj^iAtZj, / = ...N, (43) 
where Z is a Gaussian vector of N independent standard random variables. The Asian payoff is: 

1 N 

= {h{Z) - Ks)+ with /2(Z) = -ES;(Z). (44) 

;=i 

As done in the BS setting, we find the LT-based convenient directions of stratification applying 
the LT technique to the differentiable part of the payoff function of an Asian option (in this 
dynamics we only consider options on a single asset). This is done also for the other barrier-style 
Asian options, so that we have the same directions of stratification for the three types derivative 
contracts. Applying the LT decomposition the Euler scheme becomes 

N 



Sj = Sy_i +i)c{}i- Sy_i) Af + a J Sj_iAt E ^jm^m, j = 1,...N, (45) 



m=l 



the computation of the first direction of LT decomposition consists in the following steps: 
1. Compute the partial derivatives j = 1, . . . , N: 



a lAt ~ 



1 - aAi + - J -— E Ajmem 

V m=l 



/ I 



Now denote pf^ = = (l - aAt + |,y^E™=i A;„,e™) and fif\ = a^AtS~m, 

we have 

r^ = r^41,+f^,A,r. (47) 

Remark 2. The third term in a'^^ is zero, nevertheless we show its expression because the results 
below still hold when we compute the vector a.^'^ of parameters in the l-th step, where we consider 
ei = (ia^,0,...,0),/ = l,...,N. 

1—1 times 

Proposition 3. The solution of the recurrence equation ( [47I is a linear combination of the rows of A: 

vf^ = E ' 0')^'"i = 1. ■ ■ ■ . (48) 

m=l 



where the components of vector w^^' (;'), that depends on j, are: 



'{i)=f^-^Y\<'- (49) 



i=m 



The superscripts indicate the number of the direction under consideration and the proof can 
be obtained by iteration. 

Remark 3. Note that w^p {]) = jSjj'j with the assumption that Y\i£(Zi^t^ ~ ^ '^^^ ^mHj + 1) = 



2. Denote h{e) = h{Z) = /i(Ae) then 



dm 1 ^ (1) , , 



Corollary 1. ^ in equation i fjof is a linear combination of the rows of A: 

Lpf^ = Lf^ji' VNeN, (51) 



where 



As for Proposition |3| the proof can be obtained by iteration. 
Remark 4. t'^^ = ^'^p_^ = w^^ [N). 

3. The first optimal direction is established by the following theorem. 

Theorem 1. The first column of the matrix A, solution of the LT optimization problem, in the case 
of Asian options assuming the Euler discretization of the CIR model is: 

t{i) 

= W\\' ^''^ 

with t being the vector defined in Corollary^ 

Proof. Knowing that the scalar product t'^' ■ A.i attains the maximum when the two vectors 
are parallel, we can conclude that the optimal A*^ is proportional to t^^^. After normalization 
the optimum solution is given by equation 153) . □ 

Remark 5. We observe that, ifZ = 0, after some algebra, the Euler discretization is simply 

Sj-p= {1- aAt) (Sj^i - }i) (54) 

then 

Sj = (1 - izAty (So - + /' (55) 

We use the results of this remark to simplify the computational cost to find the first direction of 
stratification. 

4. In order to compute the remaining optimal columns we need to repeat the procedure illus- 
trated in steps 1 to step 3. As far as the calculation of the l-th column is concerned, one 

needs to evaluate ^^^^J^'^ and accordingly the quantities p^, cc'p , (5^ , V/, and the compo- 
nents of the vectors w''' and t*-'). All the results in Proposition [3} Corollary [1] and Theorem 
[1] remain valid while now considering the quantities with superscripts /. The orthogonal 
directions LT are then obtained by orthogonalization. 
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5.3 Linear Approximation in the Black-Scholes Market 



Hereafter we describe how to find the directions of the LA technique in the case of a BS dynam- 
ics. Since the payoff function is not differentiable, as for the LT method we consider only the 
differentiable part g — Kg. The gradient has components: 



g^g) MN r MN ' 

E '^km exp < Ffc + E C^e/ 



den 



k=l 



1=1 



then. 



V^(0) 



and in general Vg(ffl) = 



(56) 



In the above derivation we assume that C = C'^H since we do not need to introduce any orthog- 
onal matrix and the Cholesky decomposition of the autocorrelation matrix of a Brownian motion 
is explicitly known. It turns out that the LT and the LA methods return the same first order 
direction. Nevertheless, the latter approach can produce different directions changing the value 
at which the gradient is calculated. In contrast, the LT procedure admits solution only assuming 
the starting points strategy described above. Hence, the LA is more flexible and in particular 
the new algorithm does not require an incremental QR decomposition to find the new directions. 
Indeed, if we would look for orthogonal directions a unique orthogonalization would be required; 
consequently, the LA computational cost is much lower. Moreover, the mathematical derivation 
is simpler. 



5.4 Linear Approximation in the CIR Market 

We now illustrate how to apply the new LT approach for the derivative contracts above in CIR 
d5mamics. Consider the Euler discretization scheme in equation and compute the following 
partial derivatives for j,l = 1, . . . , N: 



dSj 



then 



and the gradient is 



dSjiO) 



dZi 



as,-_i(o) / 

(1 - aAt) + cr^AtSj^,{0)3ji, 



VSy(O) 



{l-ocAty~'a^/AtS^ 
{l-aAty~^cr.yAtSi{0) 



(1 - a At) a J At S 1^2(0) 



(7^AfSy_i(0) 




Due to Proposition m the LA first optimal direction is given by the normalized sum of VS^(O), / = 
1, . . . , N. Further directions are obtained by iterating this procedure with a starting points rule. 
Alternatively, we can choose the evaluation points as in the LT strategy or the components of the 
l-th direction for the starting point of the gradient for the 1 + 1-direction. 



(57) 



(58) 
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6 Numerical Illustrations 



We now illustrate the results developed in the previous sections through examples and numerical 
experiments. As mentioned before, we consider the BS and the CIR dynamics and different exotic 
path-dependent options. All the numerical procedures have been implemented in MATLAB on a 
computer with Intel Pentium M, 1.60 GHz, 1 GB RAM. In the numerical illustrations we consider 
K = 1000 strata and Ng = 2 x 10^ total number of scenarios so that for orthogonal directions 
we have a constant allocation rule (which, in this case, coincides the proportional rule as the 
strata are equiprobable) with 2000 random draws in each stratum {const in the tables). When we 
consider non-orthogonal directions the constant allocation rule is not proportional anymore since 
the strata are not equiprobable. For the optimal allocation rule {opt), the standard deviations have 
been computed by a first pilot run and then they have been used in a second stage to determine 
the stratified estimator. 

We report the estimated variances and the total computational times with constant and opti- 
mal allocation. We compare the variances employing the directions of stratification returned by 
GHS (see Glasserman et al. [8J), LT, LA, the PGA and their combination. Note that the GHS 
procedure requires the calculation of an importance sampling direction that is a computation- 
ally demanding task. In our experiments we report the variances due to the stratification only 
in order to compare the relative efficiency of the pure stratification methods. As far as the PGA 
directions are concerned, they consist of the eigenvectors associated to the highest eigenvalues of 
the autocorrelation matrix of the multi-dimensional Brownian motion that drives the BS dynam- 
ics. In contrast, since the GIR d5mamics is not Gaussian, in a first pilot run with a 2000-sample 
we compute the MG estimation of the autocorrelation matrix of the price dynamics and then 
calculate its eigenvectors and values. We employ a Euler scheme that always takes the positive 
value of the square-root term because it was shown that this exhibits the smallest discretization 
bias among Euler GIR-discretizations (see Andersen |2|). Even if this d5mamics is not normal, 
the z-th step price, given the i — 1-th one, is normal and this can justify the use of the PGA in 
the GIR dynamics. We consider the multiple combination of two directions of stratification. Our 
algorithm and considerations are also applicable to additional directions but, due to the so called 
curse of dimensionality, this would require a higher number of strata and hence a higher number 
of total samples that would considerably increase the computational burden. Finally, we compare 
these stratified estimators to LHS-based estimators (see Owen j ij] or Stein [ jiy j for more on this 
topic). Stein [17] proved that LHS eliminates the variance of the additive part of the integrand 
(payoff) function and hence produces an important variance reduction when coupled with LA 
or LT. Unfortunately, it is difficult to numerically compute the asymptotic variance in the central 
limit theorem for the LHS estimator. LHS is characterized by a fixed multiple allocation rule that 
has a high computational cost. Our purpose is to compare this very high-dimensional allocation 
rule to one with a lower dimension where we can adopt optimal allocation. In addition, the ex- 
pectation of interest E[)p(Z)] is equal to E[!/7(OZ)] where O is a general orthogonal matrix. In a 
standard MG simulation the variance of the two estimators does not depend on O but in contrast, 
the accuracy of LHS-based estimators critically depend on the choice of O. Our simulations adopt 
the orthogonal matrix produced by the LT decomposition that has been shown to be an efficient 
choice (see Sabino |i6J). 



6.1 Asian Options in the Black-Scholes Market 

Our first example is the pricing of arithmetic Asian options on a single risky security defined by 
equation |j33j with M = 1 . For simplicity we assume that the time grid is regular with time steps 
tj = iAt, i = 1, . . . ,N. This permits a simple derivation of the PGA and the Gholesky decomposi- 
tion of the autocorrelation matrix of the Brownian motion (see Akesson and Lehoczky [1]). Table 



i(a) reports the input parameters used in the simulation with different mone5mess of the options. 
We remind that in this setting LT and LA provide the same first order direction. Tables |5][7| report 
the numerical results obtained and the total computational times: all the procedures return unbi- 
ased estimates of the option prices while giving remarkably different variances. All the stratified 
techniques give a variance reduction that is particularly remarkable with the GHS and the LA (LT) 
methods. The PGA orthogonal directions (one dimensional and two dimensional) give a modest 
effect also taking into account the computational times. The main observation is that GHS and LA 
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Table i: Input Parameters in the BS d3mamics 
(a) Arithmetic Asian Options (b) Arithmetic Asian Barrier Options 



So 


Ks 


N 


r 


u 


T 


50 


45,50,55 


64 


0.05 


0.3 


1 



So 


Ks 


B 


N 


r 


a 


T 


50 


50,55 


60,70,80 


16 


0.05 


0.1 


1 



Table 2: Angles between the Stratifying Directions in degrees 



(a) Arithmetic Asian Options 





Xs =45 


Ks =50 


Xs = 55 


LA-GHS 


1.35 


1.04 


1.74 


LA-PCA 


54.62 


52.73 


51.60 


GHS-PCA 


56.60 


53.83 


53.30 



(b) Arithmetic Asian Barrier Options 





Ks-- 


= 50 


Ks-~ 


= 55 




B = 60 


B = 70 


B = 70 


B = 80 


LA-GHS 


0.37 


0.37 


0.75 


0.75 


LA-PCA 


51.95 


51.95 


51.89 


51.89 


GHS-PCA 


51.67 


51.67 


51.10 


51.10 



(LT) show the same computational cost and the same variance reduction. Both LA and GHS give 
a remarkable variance reduction, of a factor of more than 100 in the case of constant allocation and 
of several hundreds in the case of optimal allocation. However, given the parameters in Table [i(a)} 
we stress the fact that the computational time required for the calculation of the direction is really 
a small part of the total time requested for all the proposed procedures. In contrast, with a really 
high problem dimension (i.e. a dimension 1000 typical in financial applications), the solution of 
the GHS optimization problem becomes a hard task depending on the starting guess and its com- 
putational burden has a relevant influence. In contrast, the LA (LT) algorithm consists in a simple 
vector 0{N) calculation that is feasible even in high-dimensional problems. Table 2(a)| reports 
the angles (in degrees) between the discussed directions. The GHS and LA directions are almost 
parallel meaning that the GHS algorithm is not so sensitive to the moneyness and this justifies 
the equal performance in terms of variance reduction of the LA method. As mentioned before, 
the PCA direction does not furnish a relevant variance direction and hence the non-orthogonal 
2-dimensional stratification that employs such a direction always returns a lower accuracy than 
the GHS or LA methods. Moreover, the orthogonal GHS or LA bi-dimensional stratifications give 
variance reductions that are about 4 times lower than the corresponding one-dimensional ones. 
We remind that the two settings have the same number of strata so that we can conclude that the 
second order direction has a lower impact on the variance reduction and, with these directions of 
stratification, it is more efficient to employ a stratified MC estimator with a single direction. We 
conclude the study for the simple Asian options with the comparison between the accuracies of 
the LHS and the stratified sampling with a single direction with optimal allocation. The results 
shown in Table [5] illustrate that the LHS never outperforms the optimal allocation. Indeed, the 
LHS-based variance is at least two times the variance obtained with the stratified estimator with 
optimal allocation. Moreover, the computational cost is a lot higher, almost twice as high as the 
times needed for the optimal allocation. All these arguments strongly favor the use of convenient 
directions with optimal allocation. 

We modify the Asian option example by adding a knock-out barrier at expiration T or at 
each sampling date so that the option pays nothing if the asset price is above the barrier. Due to 
the discontinuous payoff of barrier options, the GHS optimization problem is a demanding task 
especially when the barrier is at each time step (indeed Glasserman et al. [8] did not elaborate this 
possibility). In contrast, the LA (LT) focuses only the continuous part of the payoff function. Table 
|i(b)| reports the input parameters used in the simulation with different moneyness and barriers. 
The values of the barriers should be larger than the strike prices but not too high otherwise the 
pricing problem would almost boil down into the case without barrier. Also for barrier options 
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Table 3: Input Parameters and Angles between Directions of Stratification for Basket Options. 

(a) Input Parameters. (b) Angles in degrees 



M 


So 


P 


r 


u 


T 


40 


Linear 20-60 


0.5 


0.05 


Linear 0.1 - 0.4 


1 





Ks =30 


Ks =40 


Ks =50 


LA-GHS 


2.76 


3.11 


2.52 


LA-PCA 


64.74 


65.04 


65.19 


GHS-PCA 


62.29 


62.02 


62.47 



(barrier at expiry), we notice that GHS and LA give directions of stratification that are almost 
parallel as illustrated in Table 2(b) This justifies the approximation of the LA method and its use 
for stratified MC to price the two types of barrier options. In addition, the GHS algorithm is not 
applicable to Asian options with a complete barrier. Different approaches should be employed 
in order to improve the stratification efficiency for barrier-style options, as suggested in Etore et 
al. |4|, but these are nevertheless computationally expensive and use orthogonal directions. The 
stratified MC does not return variances as low as for plain Asian options, especially when the 
barrier is close to the strike price. For example, the case of Asian options with barrier B = 80 
(both at expiry and at each sampling date) and with strike Kg = 55 displays a variance reduction 
of several hundreds with a computational time that ranges between 22% and 55% higher than the 
standard MC. However, when the barrier and the strike price are Kg = 50 and B = 60, respectively, 
the variance reduction is lower with an extra effort ranging from 22% and 50% with respect to the 
standard MC. 

The numerical simulation of the prices of Asian basket options with a barrier close to the 
strike price, both at expiry and at all the monitoring times, shows that stratifying along multiple 
directions can be worthwhile. Indeed, if Xg = 50 and B = 60, the multiple stratification enhances 
the accuracy of the estimation compared to the use of a single direction. In particular, the highest 
variance reduction is achieved with the choice of non-orthogonal directions (LA-PCA) with opti- 
mal allocation. In this setting the variance reduction is of an order 100, with barrier at expiry, or 
40, with barrier at each monitoring time, and is several times higher compared to the other setting 
of stratification. 

Finally, even for Asian barrier options the LHS never outperforms the technique that displays 
the smallest variance with optimal allocation. These considerations suggest that the use of mul- 
tiple non-orthogonal directions can be worthwhile. However, finding many different multiple 
directions is not a simple task. 



6.2 Basket Options in the Black-Scholes Market 

In this example the stratification estimator once more improves the accuracy of the standard MC 
method. Indeed, in the BS market, the financial features of basket options are almost the same as 
those of arithmetic Asian options. The main difference between the two is that for Asian options 
the Gaussian variables are correlated by the autocovariance matrix of a single Brownian motion 
while for basket options the dependence is measured by the covariance matrix among the asset 
returns. In addition, both payoffs contain a (weighted) average of the exponential of a Gaussian 
random vector. Table [8] shows that for all the considered exercise prices, the stratification using the 
LA (LT) with and without optimal allocation has a remarkable variance reduction comparable to 
the one given by the GHS algorithm with the same computational considerations as in the Asian 
option example. Indeed, these two directions are almost parallel (see Table [3(b)| |. The PCA-based 
direction has again a modest effect in terms of variance reduction and the stratification over a sin- 
gle linear projection produces a better accuracy than the one that exploits two directions. Finally, 
the LHS estimator neither achieves a higher variance reduction than the stratified estimator with 
a single LA direction with optimal allocation nor does it require a lower computational effort. 



6.3 Asian Options in the CIR Market 

As a last example we consider arithmetic Asian options on a single asset in a CIR djmamics whose 
depicted parameters (in Table 4(a)) are chosen in order to ensure positive prices (2af( > u^). In this 
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Table 4: Input Parameters and Angles between Directions of Stratification in the CIR dynamics. 



(a) Input Parameters 



So 


JV 


r 


a 


/' 


a 


T 


100 


64 


0.05 


1.5 


1 


0.8 


1 



(b) Angles in degrees for Asian Options 





Ks =90 


Ks = 100 


Ks = 110 


LA-LT 


1.00 


1.00 


1.00 


LA-PCA 


43.72 


43.72 


43.72 


LT-PCA 


44.24 


41.52 


41.53 



setting the LA method and the LT decomposition do not provide the same stratification direction 
and the GHS algorithm is really difficult to apply. However, as illustrated in Table [4] the directions 
returned by the LT and LA are almost parallel. In any case the derivation of the LA solution and 
its implementation are much easier Since the CIR model is neither a Gaussian nor a lognormal 
process, the PCA decomposition is not applicable. However, in order to obtain a further direction 
we estimate a PCA-like direction as explained at the beginning of this section. Tables [ojfTTI show 
that both the LA and LT algorithms give remarkable variance reductions. The best accuracies 
are obtained with the stratification along a single direction which attains a reduction of an order 
of several hundreds, both with a constant and optimal allocation rule. The extra cost for the 
computational time is only 20%. As in the BS setting, the PCA approach is less efficient and 
requires a higher computational cost due to calculation of the sampled autocovariance matrix of 
the price process. Also in this situation the solution employing two orthogonal or non-orthogonal 
directions provides a variance reduction. Unfortunately, this choice never provides an accuracy as 
precise as the one obtained by a single direction. Moreover, the use of the fixed LHS-allocation rule 
never enhances the accuracy of the simulation more than the best low-dimensional stratification 
method with optimal allocation. 

As in the BS example, we add a knock-out barrier at expiry or at each monitoring time. For 
this latter option we must chose a barrier level that is much higher than the strike price. Indeed, 
due to the high variability of the CIR djmamics, with a low barrier value the option would easily 
knock-out producing a zero-valued price. 

As already mentioned, in the example of barrier options we adopt the same convenient direc- 
tions of stratification that we would consider without the barrier since the LA and LT approaches 
do not take into account the non-differentiable part of the payoff. Tables [H and [7] illustrate the 
results of this numerical investigation. The variance reduction is not as efficient as the case 
without barrier but in contrast, the use of multiple directions improves the efficiency of the 
simulation without highly influencing the computational cost. In addition, the combination of 
non-orthogonal directions can achieve a better variance reduction. Indeed, the combination of 
LA-PCA directions (LT and LA are almost parallel) returns a variance that ranges from 10 to 30 
times lower than that with standard Monte Carlo. Moreover, this estimated variance is always at 
least equal, for Ks = B = 170 with barrier at each monitoring time, or lower than the variance 
obtained with different combinations of stratifying directions and barrier levels. 

Finally, as in all examples, the LHS sampling coupled with LT does not provide a convenient 
alternative to stratification over few directions with optimal allocation. 

7 Concluding Remarks and Future Perspectives 

In this paper we have investigated the use of convenient multidimensional directions of stratifica- 
tion in order to enhance the accuracy of Monte Carlo methods. We have discussed directions of 
stratification that are easy to derive and display variance reductions that are comparable to those 
introduced by Glasserman et al. [8J. These solutions do not require a complex calculation and 
can be applied in really high-dimensional problems without an extra cost. In contrast, the use of 
the Glasserman et al. f8] or Etore et al. [4] methods risk to be computationally unfeasible and are 
based only on orthogonal directions. Indeed, the LT and the LA directions are computed under 
convenient approximations that lead to simple matrix operations and vector norms. Moreover, 
we have proved an algorithm that allows to correctly generate Gaussian vectors stratified along 
non-orthogonal directions. Our numerical experiments demonstrate that the proposed convenient 



17 



directions return remarkable variance reductions both in BS, where the proposed techniques dis- 
play the same variance reduction as those given by GHS, and in the CIR d5mamics. In particular, 
the use of multiple non-orthogonal directions can be worthwhile for barrier style options. More- 
over, in this work we show that the use of a few convenient directions of stratification with 
optimal allocation always outperform LHS (even in its LT-enhanced form) especially in terms of 
computational burden. A natural extension would be the combination with importance sampling 
procedures like the Robust Adaptive Technique recently proposed by Jourdain and Lelong [lo] for 
Gaussian random vectors. In addition, due to its simple derivation and its affinity with the Fox's 
greedy rule (see Fox [6|), it would be interesting to investigate how to apply the LA procedure to 
derive a Quasi-Monte Carlo version of discretization schemes for stochastic volatility models like 
those proposed by Andersen {2^ and Jourdain and Sbai [11] . 
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Table 5: Results for Arithmetic Asian Options in the BS d5mamics. 





Price 






1 Dir 


2 dirs 










MC 


GHS 


LA 


PCA 


GHS 


LA 


PCA 


GHS-PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks=45 


7.02 


var 
time 


55.89 
1 


0.32 

xl.41 


0.06 
xl.51 


0.31 
xl.41 


0.06 
xl.51 


15.46 
xl.41 


11.4 
xl.51 


1.74 

xl.41 


0.61 
xl.51 


0.94 
xl.41 


0.16 
xl.51 


10.08 
xl.41 


8.66 
xl.51 


8.12 
xl.58 


0.21 

xl.68 


8.32 
xl.58 


0.19 
xl.68 


0.06 
x3.6 


Ks = 50 


4.02 


var 
time 


36.966 
1 


0.28 
xl.41 


0.04 
xl.51 


0.31 
xl.41 


0.05 
xl.51 


20.94 
xl.41 


16.18 
xl.51 


0.95 
xl.41 


0.2 
xl.51 


0.94 
xl.41 


0.12 
xl.51 


7.77 
xl.41 


6.18 
xl.51 


9.47 
xl.58 


0.21 
xl.68 


9.21 
xl.58 


0.2 

xl.68 


0.06 
X3.24 


Ks = 55 


2.06 


var 
time 


20.357 
1 


0.3 
xl.41 


0.02 
xl.51 


0.31 
xl.41 


0.03 
xl.51 


10.52 
xl.41 


7.75 
xl.51 


1.06 
xl.41 


0.28 
xl.51 


0.93 
xl.41 


0.09 
xl.51 


7.54 
xl.41 


3.8641 
xl.51 


7.4 
xl.58 


0.13 
xl.68 


7.49 
xl.58 


0.13 
xl.68 


0.06 
X3.77 



Table 6: Results for Arithmetic Asian Options with a Barrier at Expiry in the BS djmamics. 





Price 






1 Dir 


2 dirs 










MC 


GHS 


LA 


PCA 


GHS 


LA 


PCA 


GHS-PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks =50 
B = 60 


1.38 


var 
time 


2.99 

1 


1.13 
xl.41 


0.3 
xl.41 


1.13 
xl.41 


0.31 
xl.35 


2.99 
xl.41 


2.99 
xl.35 


0.54 
xl.47 


0.23 
xl.47 


0.83 
xl.47 


0.19 
xl.40 


1.24 
xl.47 


0.93 
xl.40 


0.33 
xl.50 


0.02 
xl.50 


0.32 
xl.55 


0.02 
xl.50 


1.02 
X3.91 


Ks =50 
B = 70 


1.9 


\'ar 
time 


4.8 

1 


0.13 
xl.41 


0.01 
xl.41 


0.13 
xl.41 


0.01 
xl.35 


4.77 
xl.41 


4.77 
xl.35 


0.3 
xl.47 


0.16 
xl.47 


0.15 
xl.47 


0.02 
xl.40 


1.28 
xl.47 


0.99 
xl.40 


0.41 
xl.55 


0.02 
xl.50 


0.68 
xl.55 


0.02 

xl.50 


0.13 

X3.90 


Ks =55 
B = 70 


0.19 


var 
time 


0.49 
1 


0.04 
xl.41 


0.00074 
xl.41 


0.04 
xl.41 


0.00082 
xl.35 


0.48 
xl.41 


0.48 
xl.35 


0.04 
xl.47 


0.0035 
xl.47 


0.06 
xl.47 


0.0039 
xl.40 


0.22 
xl.47 


0.06 
xl.40 


0.17 
xl.55 


0.0038 
xl.50 


0.16 
xl.55 


0.0037 
xl.50 


0.04 
X3.89 


Ks =55 
B = 80 


0.2 


var 
time 


0.55 
1 


0.0016 
xl.41 


0.00026 
xl.41 


0.0018 
xl.41 


0.00058 
xl.35 


0.55 
xl.41 


0.54 
xl.35 


0.05 
xl.47 


0.0037 
xl.47 


0.06 
xl.47 


0.0038 
xl.40 


0.22 
xl.47 


0.06 
xl.40 


0.18 
xl.50 


0.0048 
xl.55 


017 
xl.50 


0.0048 
xl.55 


0.0018 
X3.91 



Table 7: Results for Arithmetic Asian Options with a Complete Barrier in the BS d5mamics. 





Price 






I Dir 


2 dirs 










MC 


LA 


PCA 


LA 


PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks =50 
B = 60 


1.22 


var 
time 


2.42 
1 


0.85 
Xl.54 


0.23 
xl.14 


2.42 
xl.54 


2.39 
xl.14 


0.54 

xl.54 


0.12 
xl.14 


1.23 
xl.54 


0.92 
xl.14 


0.53 

xl.56 


0.07 
xl.22 


0.77 
X3.80 


Ks =50 
B = 70 


1.89 


var 
time 


4.76 
11.17 


0.14 
xl.54 


0.0047 
xl.14 


4.75 
xl.54 


4.75 
xl.14 


0.16 
xl.54 


0.02 
xl.14 


1.29 
xl.54 


1 

xl.14 


1.52 
xl.56 


0.02 
xl.22 


0.15 
X3.81 


Ks =55 
B = 70 


0.19 


var 
time 


0.47 
1 


0.041 
xl.54 


0.00087 
xl.14 


0.47 
xl.54 


0.46 
xl.14 


0.06 
xl.54 


0.0038 
xl.14 


0.22 
xl.54 


0.06 
xl.14 


0.14 
xl.56 


0.0036 
xl.22 


0.04 
X3.85 


Kc; - 55 

B - 80 


0.2 


\'ar 
time 


0.55 
1 


0.0015 
xl.54 


0.000059 
X 1.14 


0.55 
xl.54 


0.53 
xl.14 


0.05 

xl.54 


0.0038 
xl.14 


0.22 
xl.54 


0.06 
xl.14 


0.056 
x 1.56 


0.0048 
xl.22 


0.002 
X3.83 



Table 8: Results for Basket Options in the BS dynamics. 





Price 






1 Dir 


2 dirs 










MC 


GHS 


LA 


PCA 


GHS 


LA 


PCA 


GHS-PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks = 30 


11.58 


var 
time 


61.77 
1 


0.09 
xl.48 


0.06 
xl.60 


0.1 

xl.48 


0.06 
xl.60 


31.54 
xl.48 


21.63 
xl.60 


0.93 
xl.48 


0.29 
xl.60 


0.91 
xl.48 


0.25427 
xl.60 


21.17 
xl.48 


18.34 
xl.60 


6.33 
xl.75 


0.24 
xl.79 


5.27 
xl.75 


0.25406 
xl.79 


0.06 
2.87 


Ks =40 


4.15 


var 
time 


34.88 
1 


0.07 
xl.48 


0.03 
xl.60 


0.08 
xl.48 


0.04 
xl.60 


24.91 
xl.48 


17.74 
xl.60 


0.84 
xl.48 


0.15 
xl.60 


0.86 
xl.48 


0.15 
xl.60 


19.1 
xl.48 


16.71 
xl.60 


3.9 
xl.75 


0.12 
xl.79 


3.69 
xl.75 


0.13214 
xl.79 


0.1 
2.81 


iCs = 50 


0.93 


var 
time 


8.92 
1 


0.04 
xl.48 


0.004 
xl.60 


0.05 
xl.48 


0.005 
xl.60 


3.92 
xl.48 


3.88 
xl.60 


0.8 
xl.48 


0.06 
xl.60 


0.81287 
xl.48 


0.06 
xl.60 


3.05 
xl.48 


2.18 
xl.60 


2.87 
xl.75 


0.04 
xl.79 


2.55 
xl.75 


0.05 
xl.79 


0.08 
2.89 



Table 9: Resiilts for Asian Options in the CIR djmamics. 





Price 






1 Dir 


2 dirs 










MC 


LT 


LA 


PCA 


LT 


PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks = 90 


15.63 


var 
time 


427.73 
1 


1.85 
xl.2 


1.09 
xl.22 


1.54 
xl.2 


0.9 
xl.22 


115.73 
xl.5 


106.85 
xl.6 


9.3 
xl.2 


2.28 
xl.22 


51.21 

xl.5 


40.61 
xl.6 


9.13 
xl.55 


4.62 
xl.55 


1.08 
x2.76 


Ks = 100 


10.6 


var 
time 


310.11 
1 


1.49 
xl.2 


0.67 
xl.22 


1.22 
xl.2 


0.54 
xl.22 


97.22 
xl.5 


69.7 
xl.6 


8.75 
xl.2 


1.73 
xl.22 


53.03 
xl.5 


25.73 
xl.6 


8.92 
xl.55 


1.66 
x 1.554 


1.02 
x2.75 


Ks = 110 


6.95 


var 
time 


212.19 
1 


1.18 
xl.2 


0.37 
xl.22 


xl.2 


0.29 
xl.22 


82.25 
xl.5 


54.28 
xl.6 


8.72 
xl.2 


1.26 
xl.22 


40.34 
xl.5 


20.69 
xl.6 


8.29 
xl.55 


2.22 
xl.55 


0.9 
X2.76 



Table 10: Results for Arithmetic Asian Options with a Barrier at Expiry in the CIR d5mamics. 





Price 






1 Dir 


2 dirs 










MC 


LT 


LA 


PCA 


LT 


PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks = 100 
B = 110 


2.63 


var 
time 


60.43 
1 


45.76 
Xl.2 


17.77 
xl.22 


45.78 
xl.2 


17.22 
xl.22 


55.81 
xl.5 


38.69 
xl.6 


26.19 
xl.2 


9.17 
xl.22 


40.61 
xl.5 


12.49 
xl.6 


20.23 
xl.55 


3.08 
xl.55 


39.46 
x2.91 


Ks = 110 
B = 120 


1.82 


var 
time 


41.64 
1 


32.64 
Xl.2 


8.1 
xl.22 


32.55 
xl.2 


7.85 
xl.22 


38.69 
xl.5 


26.43 
xl.6 


20.76 
xl.2 


5.77 
xl.22 


26.4 
xl.5 


6.27 
xl.6 


11.95 
xl.55 


1.26 
xl.55 


28.52 
x2.87 


Ks = 100 
B = 120 


3.46 


var 
time 


81.21 
1 


34.54 
xl.2 


20.77 
xl.22 


31.61 
xl.2 


20.3 
xl.22 


53.62 
xl.5 


33.82 
xl.6 


37.05 
xl.2 


15.01 
xl.22 


50.05 
xl.5 


15.57 
xl.6 


21.19 
xl.55 


4.5 
xl.55 


48.97 
X2.87 



Table 11: Results for Arithmetic Asian Options with a Complete Barrier in the CIR d5TT.amics. 





Price 






1 Dir 


2 dirs 










MC 


LT 


LA 


PCA 


LT 


PCA 


LA-PCA 


LHS 










const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 


const 


opt 




Ks = 100 

B = 180 


2.84 


var 
time 


42.98 
1 


25.39 
xl.21 


7.44 
xl.l 


25.38 
xl.21 


7.32 
xl.l 


37.52 
xl.33 


27.98 
xl.23 


22.4 
xl.21 


6.25 
xl.l 


30.14 
xl.33 


16.5 
xl.23 


21.37 

xl.33 


5.57 
xl.23 


22.58 
x2.82 


Ks = 110 
B = 180 


1.1 


var 
time 


14.03 
1 


9.51 
Xl.21 


2.05 
xl.l 


9.59 
xl.21 


2.05 
xl.l 


12.68 
xl.33 


8.37 
xl.23 


8.63 
xl.21 


1.73 
xl.l 


11.33 
xl.33 


4.76 
xl.23 


8.35 
xl.33 


1.58 
xl.23 


8.79 
X2.85 


Ks = 100 
B = 170 


1.79 


var 
time 


23.7 
1 


15.86 
xl.21 


4.65 
xl.l 


15.96 
xl.21 


4.58 
xl.l 


21.59 
xl.33 


15.21 
xl.23 


10.75 
xl.21 


3.44 
xl.l 


15.97 
xl.33 


8.75 
xl.23 


13.73 
xl.33 


3.47 
xl.23 


15.08 
X2.79 
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